Variational theory of two-fluid hydrodynamic modes at unitarity 
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We present the results of a variational calculation of the frequencies of the low-lying Landau 
two-fluid hydrodynamic modes in a trapped Fermi superfluid gas at unitarity. Landau's two-fluid 
hydrodynamics is expected to be the correct theory of Fermi superfluids at finite temperatures close 
to unitarity, where strong interactions give rise to collisional hydrodynamics. Two-fluid hydrody- 
namics predicts the existence of in-phase modes in which the superfluid and normal fluid components 
oscillate together, as well as out-of-phase modes where the two components move against each other. 
We prove that at unitarity, the dipole and breathing in-phase modes are locally isentropic. Their 
frequencies are independent of temperature and are the same above and below the superfluid tran- 
sition. The out-of-phase modes, in contrast, are strongly dependent on temperature and hence, can 
be used to test the thermodynamic properties and superfluid density of a Fermi gas at unitarity. 
We give numerical results for the frequencies of these new modes as function of temperature in an 
isotropic trap at unitarity. 

PACS numbers: 03.75.Kk, 03.75.Ss, 67.40.-w 
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I. INTRODUCTION 

Landau's two-fluid hydrodynamics [l], 0] is the theory 
of the finite temperature dynamics of all superfluids (with 
a two-component order parameter) when collisions are 
sufficiently strong to produce a state of local thermody- 
namic equilibrium. Recent experiments have begun to 
probe the collective modes in trapped superfluid Fermi 
gases with a Feshbach resonance H, Q}. At unitarity, 
the magnitude of the s-wave scattering length a s that 
characterizes the interactions between fermions in differ- 
ent hyperfine states diverges (|a s | — > oo). Owing to the 
strong interaction close to unitarity, we expect that the 
dynamics of superfluid Fermi gases with a Feshbach res- 
onance at finite temperatures are described by Landau's 
two- fluid hydrodynamic equations 

Solving Landau's two-fluid equations for trapped gases 
is difficult due to the fact that the density profiles of 
the superfluid and normal fluid components are highly 
nonuniform, making a reliable "brute-force" numerical 
calculation very challenging [1, 0] • In a recent paper [j| , 
an alternate variational formulation of Landau's two-fluid 
equations was developed. Following the approach pio- 
neered by Zaremba et al. Q, we use a simple ansatz for 
the superfluid and normal fluid velocity fields based on 
exact solutions at T = and above T c . This gives alge- 
braic equations for the variational parameters describing 
the breathing and dipole two-fluid modes. The coeffi- 
cients in these equations involve spatial integrals over 
equilibrium thermodynamic quantities. This approach is 
simpler than solving the two-fluid equations directly for 
trapped gases. In the present paper, we report numerical 
results for the breathing and dipole mode frequencies at 
unitarity for an isotropic trap based on this variational 
method. However, our general approach can also be used 



away from unitarity. 

We discuss the in-phase breathing mode at unitar- 
ity since this mode has been studied extensively in re- 
cent experiments @, H E|. In particular, wc examine 
the surprising results of the experiments by Thomas and 
coworkers Q that have shown the frequency of this in- 
phase mode to be almost independent of temperature, 
remaining within a few percent of its T = value even 
well above the superfluid transition temperature T c . Our 
analysis of the Landau two-fluid equations at unitarity 
shows that the in-phase breathing and dipole hydrody- 
namic modes are locally isentropic, mode, with the su- 
perfluid and normal fluid moving with the same velocity, 
v s (r, t) = v n (r, f). We find that the frequencies of these 
in-phase modes are independent of temperature, given by 
their T = value at all temperatures. 

Of greater interest are the out-of-phase breathing and 
dipole modes, which have not been studied experimen- 
tally. These modes involve an oscillation of the trapped 
superfluid where the superfluid and normal fluid compo- 
nents move against each other, in contrast to the in-phase 
modes where these components move together. The out- 
of-phase modes are predicted to be strongly temperature- 
dependent and should provide a useful tool to test the 
microscopic model used for the thermodynamic proper- 
ties. 

In a companion paper [io| . we show how these two- 
fluid modes can be measured using standard two-photon 
Bragg scattering techniques |Tl[. Extending the varia- 
tional method described in this paper, we show the den- 
sity response function has resonances at the breathing 
and dipole mode frequencies. 

In our variational theory [f| , calculation of the frequen- 
cies of the two-fluid modes requires knowing the values 
of a number of thermodynamic quantities. At unitar- 
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ity, however, the variational equations simplify with only 
two thermodynamic quantities required for the dipole 
and breathing mode frequencies: the superfluid density 
p s and the isentropic compressibility (dp/dp) s . In this 
paper, we calculate the latter quantity at unitarity us- 
ing the fluctuation theory developed in Ref. [12j, which 
is an improved version of the original theory of Nozieres 
and Schmitt-Rink (NSR) As shown in Ref. [3, this 
theory gives thermodynamic quantities at finite tempera- 
tures which are in excellent agreement with ab-initio cal- 
culations [3, [TBI, [l6| and recent experimental measure- 
ments [TtJ • The superfluid density we use is also based 
on the NSR fluctuation theory Ojj]. The spatially- 
varying compressibility and superfluid density that enter 
our variational two-fluid equations are calculated within 
a local density approximation (LDA) using our results 
for a uniform Fermi superfluid. 

He et al. Q have also reported results for the two-fluid 
modes in an isotropic trap, based on a direct numerical 
solution of the Landau two-fluid differential equations. 
While there is some ambiguity in identifying the nature 
of the oscillations in Ref. [7| , the in-phase breathing mode 
is found to be temperature-independent, in agreement 
with our variational results. However, the temperature 
dependence of the out-of-phase mode breathing mode is 
very different from what we obtain (see Section IVII[) . 

Heiselberg [2(| has discussed the first and second sound 
velocity in the BCS-BEC crossover for a uniform gas. 
In this case, the solutions of the two-fluid equations are 
known (plane waves). For the thermodynamic functions 
which are needed, Heiselberg worked these out in the 
BCS and BEC limits and interpolated these results to 
describe the unitarity region. Our work makes a major 
extension of this previous study since we deal with a non- 
uniform trapped superfluid and use a microscopic theory 
for the thermodynamic functions and the superfluid den- 
sity at unitarity. 

In Section |TT1 we discuss some of the features of "uni- 
versal" thermodynamics valid at unitarity [2l[. We use 
these results in Section UIT1 to prove that Landau's two- 
fluid hydrodynamic equations predict a locally isentropic 
breathing mode at unitarity, corresponding to a situa- 
tion where both the normal and superfluid components 
move with the same local velocity. In Section IIV1 we 
review the variational formulation of Landau's two-fluid 
equations given in Ref. [5]. In Section El we discuss 
the NSR results for the temperature dependent isentropic 
compressibility and superfluid density which we need as 
inputs in our variational solutions. In Section fVTl we re- 
formulate the equations for the breathing modes derived 
in Ref. [5[ in a more useful form for use at unitarity. 
In Section IVII1 we show that the predictions of universal 
thermodynamics allow us to derive simple expressions for 
the breathing mode frequencies at unitarity. Numerical 
results for the temperature dependence of the frequency 
of the out-of-phase breathing mode are also given for a 
trapped gas using a local density approximation (LDA). 
In Section lVIIIl we calculate the temperature dependence 



of the out-of-phase dipole mode frequency. 

In Appendix fXl we compare the isentropic breathing 
mode in trapped Fermi superfluid gases with first-sound 
in superfluid 4 He, which is also a locally isentropic mode. 
Appendices [B] and [C] discuss the low and high tempera- 
ture limits of the frequency of the out-of-phase breathing 
mode using a BCS mean-field theory (without fluctua- 
tions). These calculations confirm the main features of 
the LDA results given in the text, still within the same 
variational ansatz. 



II. THERMODYNAMICS AT UNITARITY 

In this Section, we review the features of universal ther- 
modynamics at unitarity plj and use these to derive a 
number of thermodynamic identities at unitarity that will 
be used throughout this paper. 

In a dilute, uniform system of interacting fermions, 
there are three microscopic length scales (for a recent re- 
view and references on Fermi gases, Giorgini, Pitaevskii, 
and Stringari [13]). The three length scales are the mean 

interparticle spacing n^ 1 ' 3 , the thermal wavelength Xp = 
2tt /rnksT (throughout this paper we set h = 1), and 
the s-wave scattering length a s that completely charac- 
terizes the interaction between different species (denoted 
by the f, j.) of fermions in the low-density limit. Here, 
up = (2mei?) 3 / 2 /37r 2 is the density of both species of 
fermions (i.e., np = n-f + njj, where cf is the Fermi en- 
ergy of an ideal gas. The corresponding energy scales are 
the kinetic energy ep, ksT, and the interaction energy 
(which can be expressed as a functional of the density 
uf and a s ). At unitarity, the scattering length diverges, 
meaning that the only remaining length scales are the 

— 1/3 

interparticle spacing n F and the thermal wavelength, 
as first argued by Ho [2l|. This also implies that at uni- 
tarity, the only energy scales are the Fermi energy and 
ksT. Consequently, the only dimensionless energy scale 
at unitarity is ksT/ep = fcsT '/fcsTp . This immediately 
means that all thermodynamic functions at unitarity can 
be written in dimensionless form as a function of the ratio 
T/Tp. These features can be used to derive useful identi- 
ties involving the internal energy, entropy, and chemical 
potential. 

Owing to the fact that there is only one dimensionless 
energy scale, given by ksT/ksTp^p), the internal energy 
density U in a trapped Fermi gas takes the form [SJ, |2l| 

f/= pe £ (p) /s 
m 

Also, the total entropy S of a fluid element of small (in- 
finitesimal) volume AV is 0, [2l[ 

S = Nk B f s [T/Tp(p)}. (2) 

Here fp and fg are dimensionless functions of the re- 
duced temperature T/Tp(p). ejr(p) is the local Fermi 
energy and is a function of the mass density p(r). N(r) = 
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p(r)AV/m is the total number of fermions in the small 
volume AV centered at position r. We emphasize that 
both the energy density U(r) and the entropy S(r) of 
a small fluid volume centered at r depend on position 
through the Fermi energy e F {p) and the local mass den- 
sity p(r). 

The total local energy density is given by E = U + 
pVext, where 



(3) 



is the harmonic trapping potential divided by the mass. 
It is standard for Landau's two fluid equations to be given 
in terms of the mass density p = ran, instead of the 
number density n. Thus we use this scaled harmonic trap 
potential. At unitarity, the energy of the small volume 
AV is thus 

E AV = Ne F (p)f E [T/T F (p)} + NmV cxt . (4) 
The pressure P is defined by 



( (E AV) \ 
{ dAV ) 



(5) 



N,S 



From Eq. @, we see that holding N and S constant 
requires holding the reduced temperature constant as 
well 131 . Thus we find 



dAV = N \ -ttxtt ) fE[T/T F {p)} 



dAV 
p 2 dep(p) 



'N,S \ /N 

" -f E [T/T F {p)] 



m dp 
2 pe F (p) 
3 m 



f E [T/T F (p)}. (6) 



Thus, at unitarity, the pressure and energy density are 
related by @ 



P=\^f E [T,T F{ p)] = \u, 
dm 6 



(7) 



the same relation one obtains in a nonintcracting Fermi 
or Bose gas. The temperature is defined by 



T = 



= ( d JL\ 



(8) 



where s — S/ AV is the entropy density. Using this, 
Eq. ||7|) implies that 



d_ fdP\ _ 2dTo _ 



dxi \ds J 3 da 



= 0, 



(9) 



since the equilibrium temperature T is spatially uniform, 
even in a harmonically confined gas with nonuniform den- 
sity. 



The chemical potential per unit mass is given by 0] 

(10) 

Combining this expression with Eq. ([7]), we also obtain 

(11) 



dP\ 2 



d P J s - 3 [ ^ y ° xt] 



Using this, we find 

d fdp\ 2dp 2au cxt 



dxi \dp J 3 dxi 3 dxi 



-y iXi . (12) 



Here we have made use of that fact that, like the tempera- 
ture T , the equilibrium chemical potential po is spatially 
uniform, V/io = 0. 

We will make use of the identities derived in this Sec- 
tion (for a Fermi gas at unitarity) throughout this paper. 



III. LOCALLY ISENTROPIC DYNAMICS 

Before discussing our variational solutions of the two- 
fluid equations in Section IIVI we use the results of Sec- 
tion [IT] to discuss some general features of the solu- 
tions of the Landau two-fluid hydrodynamic equations 
for trapped superfluid gases. In particular, Thomas et 
al. [§] argued that the (in-phase) breathing mode at 
unitarity obeys a single Euler equation for the velocity 
v = v s = v„ on the grounds of locally isentropic hydrody- 
namics. It followed from the analysis of this Euler equa- 
tion that the frequency of the breathing mode would be 
independent of temperature. This surprising result was 
consistent with their experimental results for the breath- 
ing mode. We now derive this starting from Landau's 
two-fluid hydrodynamic equations. 

We start with the continuity and conservation of en- 
tropy equations of Landau two-fluid hydrodynamics Q , 



and 



| + V.(,v„) = 0. 



The total mass current 



j = p s v s + p n v„ 



(13) 



(14) 



(15) 



is given in terms of the superfluid and normal fluid ve- 
locities v s and v„, as well as the superfluid and normal 
fluid densities, p s and p n . The sum of the superfluid 
and normal fluid densities gives the total mass density, 
p = p s + p n - The continuity equation in Eq. (f]~3|) ex- 
presses mass conservation and is always valid. Equa- 
tion (TT4"]) assumes that the entropy of the fluid is carried 
by the normal fluid and is conserved. These equations 
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describe reversible flow without any dissipation arising 
from transport coefficients ||. 

An oscillation is locally isentropic if the entropy per 
unit mass s(r,t) = s(r,t)/p(r,t) = S(r, t)/p(r, t)AV 
does not change in time as the mass element p(r,i)AU 
moves with the fluid. Defining the Lagrangian derivative 



D d „ 

— = 1- v • V 

Dt ~ dt + ' 



(16) 



locally isentropic hydrodynamics corresponds to the sit- 
uation where 



D£ 
Dt 



= 0. 



Using Eqs. (fl3|) and (fT4| . one can show that 

ds __ s_ , , 

— + v„ • Vs = -V • p s (v s - v„). 
dt p 



(17) 



(18) 



This result confirms that the dynamics of a fluid are lo- 
cally isentropic when v s — v n = v. 

For locally isentropic fluid flow, Landau's expression 
for the current in Eq. (|15p reduces to j = (p s +p n )-v = pw . 
Using this result in Landau's equation of motion for the 
current [5| 



dt 



-VP - pWoxt - p s v s • Vv s - p n v n ■ Vv r , 



-v s V • (p s v s ) - v„V • (/0„v n ), 



(19) 



it reduces to 
dt 



-VP - pVU cxt - pw ■ Vv - vV • j. (20) 



Combining this equation with the continuity equation 
given by Eq. p3[) . we obtain the following equation of 
motion for the velocity v: 



9v 
~dt 



VP 

p 



(21) 



This is precisely Euler's equation for an ideal irrota- 
tional (such that Vv 2 = 2v • Vv) fluid [23j |. general- 
ized to include the effects of an external trapping poten- 
tial. This result shows that for the special case where 
v s (r,t) = v„(r, t), Landau's two-fluid hydrodynamic 
equations reduce to Euler's equation for an irrotational 
velocity field. 

Our present discussion shows the equation of motion 
considered in Ref. [9| is a rigorous consequence of Lan- 
dau's two-fluid equations for locally isentropic flow. We 
now derive a condition for a locally isentropic (v s = v„) 
normal mode solution of the Landau two-fluid equations 
to exist. 

The linearized continuity and entropy conservation 
equations [given by Eqs. (flU]) and fT4")) ] are 



dSp 
~dt 



+ V • {p s0 V s + PnOV„) = 



(22) 



and 



dSs 



V • (s v„) = 0. 



(23) 



Introducing the displacement fields [UHI u s ,u„, 



i ,n du s {r,t) du n (r,t) 
v s ( r > *) = 51 > v « ( r » *) = 5T — . ( 24 ) 



dt 



dt 



the linearized continuity and entropy conservation equa- 
tions can be expressed in terms of these fields as 



8p(r, t) = - V ■ [p s0 (r)u s (r, t) + p„ (r)u„(r, t)] (25) 



and 



6s(r,t) = -V- [s (r)u„(r,t)] 



(26) 



These expressions will be used in deriving the conditions 
for a locally isentropic mode to exist. 

Since each mass element evolves at constant entropy in 
a locally isentropic flow, these elements do not exchange 
heat with their surroundings and hence the temperature 
remains unchanged throughout the fluid. From the lin- 
earized Landau two-fluid equations for the superfluid and 
normal fluid densities (see Eqs. (38) and (39) in Ref. [5j]), 
one can show that 



d(v s - v„) 
dt 



PnO 



VST. 



(27) 



This implies V<5T = when v s = v„, showing that the 
temperature remains constant everywhere for a locally 
isentropic mode. Thus, a locally isentropic mode is also 
a locally isothermal mode. Using 5T = (dT/ds) p 8s + 
(dT/dp) s Sp and Eqs. ([23]) and (|2"6")> . we can write the 
condition VST = as 



-) V.^u, 



dr\ 

ds) l 



V • (s u) 



= 0, (28) 



where u s = u n = u. 

To make contact with the results of Section HH we ex- 
press Eq. (|28p in terms of derivatives of the pressure. 
The pressure can be expressed in terms of the equilib- 
rium thermodynamic identity ||, 



P=-U -pV cxt +Ts + pp. 



(29) 



Treating P, T, and p as functions of the independent vari- 
ables p and s, using the Maxwell relation 



dT\ = (dp\ 



and Eqs. 



dp) s \dsj p ' 
and (Unj), one can show that 

dp 



dP 
dp 



'dp 

= "' i[ jTp 



•so 



ds 



(30) 



(31) 
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and 



^ = pj?l 



ds J 



dp 



■so 



dT 
ds 



(32) 



The gradient of the equilibrium temperature To can be 
written as 



VT 



dT 
dp 



dT\ 
9s ) p 



v Po + ( ^ ) Vs - 0. (33) 



This gives the following useful identity for a trapped gas: 



/<9T N 






Up, 




m, 



(34) 



p J 



Using Eqs. ([32]) and 1(54}. the condition in Eq. (2S) can 
be rewritten in the useful form 



\ds J 



ds J 



0. (35) 



Equation (|35jl thus gives the condition for there to exist 
a locally isentropic (or isothermal) normal mode solution 
of the Landau two- fluid equations. This relation is com- 
pletely general for an oscillation described by u. It is not 
restricted to the case of a superfluid in a harmonic trap 
at unitarity, although this is the region of interest in this 
paper. 

At unitarity, the second term in Eq. (|35p vanishes in 
accordance with Eq. © . Thus we conclude that a locally 
isentropic mode (v s = v„) exists at unitarity if either 



dP\ _ 2 
~to) P ~ 3 



T = 0, 



or if 



V(V-u) = 



(36) 



(37) 



is satisfied. The first condition given by Eq. (]36|) is triv- 
ially satisfied at T = 0. Here the normal fluid van- 
ishes and hence all particles move with the same velocity 
v s = v, and of course any oscillation will be locally isen- 
tropic. In order for a locally isentropic mode to exist at fi- 
nite temperatures, Eq. (|37|) must be satisfied. This is sat- 
isfied by the scaling solution [24] v(r, t) cx r cos ujt [equiv- 
alently u(r, t) oc rcoswi] of the hydrodynamic equa- 
tion in Eq. (f2"Tj) that describes the breathing mode. It 
is also satisfied by the generalized Kohn mode (the in- 
phase dipole mode) that we discuss in Section IVIIH This 
suggests that the existence of a purely locally isentropic 
mode is not a universal feature of hydrodynamics at uni- 
tarity, but rather is a special feature in a harmonically 
confined gas. 

In Section IVII1 we confirm that our variational solu- 
tion of the two-fluid equations gives a locally isentropic 
breathing mode with a frequency independent of tem- 
perature. We call this breathing mode the "in-phase" 



breathing mode since the normal and superfluid compo- 
nents move together, v s = v„. This is the mode studied 
by Thomas and coworkers 0. In addition, our varia- 
tional solution also predicts an out-of-phase breathing 
mode which is not locally isentropic, with a frequency 
very strongly dependent on temperature. 

In superfluid 4 He, first sound also describes a locally 
isentropic mode, a fact accounted for by Eq. (|3"5")) . How- 
ever, first sound in uniform superfluid 4 He is locally isen- 
tropic for different reasons than the in-phase breathing 
and dipole modes in a trapped Fermi superfluid at uni- 
tarity. This is discussed in Appendix [XJ 



IV. VARIATIONAL SOLUTION OF THE 
TWO-FLUID EQUATIONS 

While the preceding analysis showed that the Lan- 
dau two-fluid equations at finite temperatures admit a 
class of analytic solutions at unitarity [corresponding to 
V(V ■ u) = 0], these solutions only describe the in-phase 
[u s = u„ = u] dipole and breathing mode oscillations. 
The out-of-phase solutions of the two-fluid equations can- 
not be obtained using such a simple analysis. We shall 
use a variational method to derive expressions for the fre- 
quencies of these out-of-phase modes. In this section, we 
review the variational formulation of Landau's two-fluid 
equations developed in Ref. 5]. 

In 1950, Zilsel [25[ introduced a phenomenological ac- 
tion S[s, p, p n , v s , v„] as a function of the entropy density 
s, the total density p = p n + p s , the normal fluid den- 
sity p n , as well as the superfluid v s and normal fluid 
v„ velocities. By construction, the variation of this ac- 
tion with respect to these variables generates the Landau 
two-fluid equations. In order to generate the linearized 
two-fluid equations (the solutions of which determine the 
spectrum of normal modes), the action is expanded in 
powers of fluctuations (5p, ds, <5v s , d~v n ) about the equi- 
librium values (Pso, so, v„o, v„o) up to quadratic order. 
We assume that v s0 = and v n0 = 0, so that Sv n = v n 
and <5v s = v s . The terms in the action that describe fluc- 
tuations Sp n in the normal fluid density can be shown to 
be higher-order || and are thus neglected. The result- 
ing action describes the hydrodynamic fluctuations. It 
is further simplified by replacing the entropy and den- 
sity fluctuations ds and Sp in terms of the superfluid and 
normal fluid velocities. This can be done using the lin- 
earized continuity and entropy conservation equations in 
Eqs. (22) and (23). 

Using Eqs. (21), (25), and (25), the action that de- 
scribes hydrodynamic fluctuations (8p, 8s, v s , v n ) can be 
expressed in terms of the two displacement fields u s and 



S" 2 



drdt I -p s0 u 2 s 



1 



• 2 

PnOU„ 



6 



dT 
dp 



[V • (s U„)] [V • (PsQU s + p n0 Un 



-i(f)[v.( S0 „jf 



(38) 



Here p is the chemical potential per unit mass defined in 
Eq. ([TO]) and T is the temperature defined in Eq. f8j. 

Formulating the linearized two-fluid equations in terms 
of the variation of an action as in Eq. ([55]) allows us 
to develop variational solutions of these equations by 
making an ansatz for the displacement fields u s (r,t) and 
u„(r, t). This was done in Ref. ||, extending earlier work 
in Ref. |8| for the two-fluid modes of a trapped Bose- 
condensed gas at finite temperatures. Our variational 
ansatz for each Cartesian component of the displacement 
fields is 



u si (r,t) 
u ni (r,t) 



asifi(r) cos ut, 
a n igi{r) cos tot. 



(39) 



The constants a S i and a n i are the variational parameters. 
With an ansatz of this form, the variational equations 
reduce to 



da. 



da r . 



(40) 



Once some suitable ansatz is made for the functions /i(r) 
and <7i(r) in Eq. (|3"9")l . these equations can be used to 
generate variational solutions of the two-fluid equations 
and the corresponding normal mode frequencies to. 

For gases confined in a harmonic trap, there exist sim- 
ple trial functions for /^(r) and gi(r) which are sufficiently 
close to the exact solutions that good results for the mode 
frequencies to are obtained by considering only a single 
expansion term as in Eq. (|3"9")) [8|. The choice of ansatz 
for the displacement fields at finite temperatures used in 
Ref. [|[ for the dipole and breathing modes are guided by 
the known exac t hy drodynamic solutions at T — [26[ 
and T > T c p7l . l28j |. For the breathing mode, we use 



/i(r) = Xi, gi(r) = x { 



(41) 



For an isotropic trap, the breathing mode in Eq. (|39D is 
described by a S i = a s and a n i = a n , in which case the 
displacement fields are given by 

u s (r, t) — a s rcosu)t, u„(r, t) = a n rcosujt. (42) 

The dipole mode is characterized by displacements of 
the centre-of-masses of the two fluids along one of the 
axes of the harmonic trap, say the z axis. In this case, 
we use the following ansatz for the displacement fields: 



/ z (r) = a s , g z {r) = a n , 



(43) 



where a s and a n describe the displacements of the centre- 
of-masses of the two fluids from the trap centre. This 
ansatz describes a uniform displacement field, 

u s (r,i) = a s zcosu>t, u„(r,i) = a n zcosuit. (44) 



At T = where the normal fluid component van- 
ishes, the ansatz used above for the breathing and dipole 
modes are exact solutions of the quantum hydrodynamic 
equations. Similarly, for T > T c , where the superfluid 
component vanishes, u n (r, t) = a n rcosu}t and u n (r,i) = 
a„zcoswi are both solutions of the collisional hydrody- 
namic equations [13, HH . 

We expect that the ansatz given above will be a good 
approximation to the exact solutions in the superfluid 
two-fluid region. We note, however, that it is straight- 
forward to improve the results presented in this paper 
by extending our variational ansatz using a generalized 
Rayleigh-Ritz expansion Q. For the breathing mode, for 
instance, this would take the form 



N 

a^ n r r 2j cos ut. 

3=0 



(45) 



In addition to improving our numerical results for the 
lowest breathing mode (n = 1,1 = 0) frequency, this 
ansatz also allows us to solve for the higher-order (n > 
1,1 — 0) "monopole" modes (see, for instance, Ref. [29]]). 

We note that the ansatz for our breathing mode in 
Eq. (gU) satisfies V(V • u) = 0. Similarly, the dipole 
mode ansatz in Eq. (|4"4")l satisfies V • u = 0. We recall 
from our analysis in Section lrTTl that the Landau two-fluid 
equations thus require the resulting in-phase breathing 
mode to be locally isentropic (corresponding to a s = a n ) 
only at unitarity, while the in-phase dipole mode is locally 
isentropic everywhere. In Section I VIII we confirm that 
our variational solution of the in-phase breathing mode is 
described by a s = a n at unitarity (and only at unitarity). 
That a s = a n for the in-phase dipole mode is always 
correct has already been shown in Refs. [1, Q. 



V. SUPERFLUID DENSITY AND ADIABATIC 
COMPRESSIBILITY OF A UNIFORM FERMI 
GAS AT UNITARITY 

In later sections, we show that our variational solutions 
for the two-fluid dipole and breathing modes at unitarity 
require as input only two thermodynamic quantities: the 
superfluid density p s and the adiabatic compressibility 
(dp/dp) s . In this section, we discuss the approximations 
used to evaluate these quantities. 

The adiabatic compressibility (dp/dp) s can be ex- 
tracted from the equation of state for a uniform Fermi 
gas at unitarity. For this purpose, we express the chem- 
ical potential and the entropy in terms of dimcnsionlcss 
functions as a function of the reduced temperature, 



£f (p) 



■in 



U[T/T F {p)\. 



and 



phi 



m 



■f a [T/T F (p)} 



(46) 



(47) 
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FIG. 1: (color online) Superfluid density fraction in a uniform 
Fermi gas at unitarity as a function of temperature. The 
different theoretical predictions are discussed in the text. 



where the dimensionless functions / M and / s , and their 
derivatives may be calculated numerically using the fluc- 
tuation theory discussed in Refs. [12. Il7j. 

Using Eq. (|46j) . the compressibility is given by 



dp 2e F (p) , e F (p) £l d[T/T F (p)] 

7T = q - ~ — Jt* "i Z~ a ' ( - 48 ' ) 

op 3 rap m * op 

where /' = df/dT' with V = T/Tp. From the expression 
in Eq. (|47p . we see that keeping the entropy constant in 
evaluating Eq. (|48|l amounts to requiring that 



d[T/T F {p)} 
dp 



111 

Pf's 



We thus obtain 

dp 
dp 



CF (p) 

mp 



~.fu 



fufs 



ft 



(49) 



(50) 



This quantity is straightforwardly evaluated using the 
values and f s obtained from the finite temperature 
equation of state of a uniform superfluid. 

The determination of the superfluid density in the 
BCS-BEC crossover is more subtle. It has been recently 
calculated for a uniform system including Gaussian NSR 
fluctuations [IH [l9j , and by Akkineni et al. using path- 
integral Monte Carlo (PIMC) simulations [30(. We sum- 
marize these results in Fig. [T] Neither the NSR fluctua- 
tion or the PIMC calculations give results that are accu- 
rate near the superfluid transition temperature T e . The 
PIMC calculation suffers from the negative-sign problem 
for fermions, and the results are thus restricted to a small 
number of total particles N = 20. The NSR-type Gaus- 
sian fluctuation theory, on the other hand, suffers from 
a re-entrance problem close to T c . This problem first 
appears around (fcj?«s) _1 = —0.5 on the BCS side and 
persists into the BEC side of unitarity . This spurious 
first-order phase transition [l9( is due to the NSR Gaus- 
sian treatment of pairing fluctuations used to calculate 



Ao and p self-consistently. The problem is equivalent 
to one that arises in a self-consistent calculation of the 
condensate density and chemical potential close to T c in 
Bose gases using the Bogoliubov-Popov approximation 
(for further discussion and references, see p. 34 of Shi 
and Griffin [3l|). However, as seen in Fig. HJ both the 
NSR and PIMC data are in good agreement at low tem- 
peratures. 

To overcome the lack of an accurate p s calculation near 
T c , we use two different sets of data for the superfluid 
density in our calculation of the out-of-phasc breathing 
and dipole modes: a fit to the NSR data from Ref. |19j 
and a scaled BCS mean-field superfluid density. A su- 
perfluid with a two-component order parameter (and a 
bosonic fluctuation spectrum) undergoes a second order 
phase transition with a superfluid density that varies as 
p s oc [T c — T) 2 / 3 close to the transition temperature, in- 
dependent of the interaction strength [32] . Our fit to the 
NSR fluctuation data thus assumes a curve of the form 
(T c - T) 2 / 3 in the region T > 0.187>. This leads to the 
fitting curve p s /p = 4.51(0.237 - T/T F ) 2 / 3 for the high 
temperature data, where T c is 0.237Tf in this fitting. 
We scale the temperature dependence of this data using 
T -> (0.225/0. 2368)T, so that the transition tempera- 
ture for a uniform gas is T c ~ 0.225Tf, as given by NSR 
theory. The final result is plotted in Fig. Q] ("scaled fit 
to NSR data"). The original NSR data from Ref. jH| is 
denoted by the blue circles (the data points used in the 
curve fitting are given by the filled circles). 

The PIMC results for p s shown in Fig. Q] have not been 
rescaled to the T c used for the other predictions. Akki- 
neni et al. {3Cj have used finite-size scaling procedures to 
obtain a T c ~ 0.25Tf- One can ignore the PIMC data 
points above 0.25TV and introduce a smooth extrapola- 
tion of the lower temperature points to vanish at 0.25Tf- 
When plotted in Fig. Q] using a rescaled T c of 0.225Tf, 
the resulting PIMC results are in fairly good agreement 
with our fitted NSR results. 

The NSR-type theories developed in Refs. [H, [U GJ, 
fioj ] includes the contributions from the BCS Fermi excita- 
tions plus the bosonic pairing fluctuations. As discussed 
in Refs. [18j, [l9(, one finds that the normal fluid density 
Pn(— P — Ps) reduces precisely to the expected Landau 
formulas in both the BCS and BEC limits. That is, the 
normal fluid is expressed in terms of Fermi excitations 
(BCS) or Bogoliubov-Popov Bose excitations (BEC), re- 
spectively. Obtaining both limits correctly is very impor- 
tant in any acceptable theory of the superfluid density in 
the BCS-BEC crossover. As noted above, however, even 
though our expression for the superfluid density reduces 
to the Landau expression on the BEC side of unitarity, 
our results are still unreliable near T c because the spec- 
trum of Bogoliubov-Popov excitations that determines 
p s is evaluated using the values of A (T) and p(T) de- 
termined self-consistently in the Gaussian NSR theory. 
The self-consistent determination of Ao(T) and p(T) in 
this approximation is equivalent to the calculation of the 
condensate density n c (T) and chemical potential p(T) 
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for a Bose gas with a Bogoliubov-Popov excitation spec- 
trum [l8|]. It is well-known (see Shi and Griffin (3l|) that 
the latter problem predicts a spurious first-order phase 
transition because one is trying to determine the conden- 
sate depletion self-consistently from the thermal excita- 
tion of collective modes with a spectrum that depends on 
the condensate fraction. 

As we have noted, the NSR-type treatment of fluctu- 
ations appears to give excellent results for the thermo- 
dynamic functions in the BCS-BEC crossover when one 
compares them with ab-initio calculations. The NSR the- 
ory does have a problem near T c near unitarity and on 
the BEC side of the crossover as a result of only consid- 
ering Gaussian fluctuations. However, we consider it the 
best available theory for the superfluid density p s in the 
BCS-BEC crossover at the present time. 

In addition to the fitted NSR data for p s , we also use a 
scaled mean-field BCS superfluid density fraction. This 
is obtained by a linear compression of the horizontal axis 
of the BCS superfluid density, 



^scaled 



[T]=P. 



rpNSR 

BCSt±c T i 



(51) 



Here Tf SR ~ 0.2257> and T^ cs ~ 0.497T F arc the 
transition temperatures of the uniform Fermi gas given 
by the NSR fluctuation theory [l2|, [l9[ and the mean- 
field BCS theory, respectively. This data is shown in 
Fig. [I] by a dotted line. While calculations [19( show 
that the BCS result for p s is only a good description 
of the superfluid density on the BCS side of resonance 
when (kpag)^ 1 < —0.5, much of the effect of "beyond 
mean-field fluctuations" is included by the scaling of T c 
in Eq. (|5ip. It is well-known that near T c a mean-field 
BCS type theory has a p s oc (T c — T) outside the region 
where fluctuations are important. 

He, Chien, Chen, and Levin [33| have also recently cal- 
culated the superfluid density in the BCS-BEC crossover 
using the pseudogap theory [34| . For comparison, in 
Fig. [JJ we also plot the pseudogap result for p s at uni- 
tarity. This is obtained by evaluating the expression 
given in Eq. 70 of Ref. [33[ using values of the gap 
and chemical potential obtained by solving Eqs. (34)- 
(36) of Ref. [23|. The pseudogap expression for p s 
is given by p s = A 2 sc /A 2 pf cs (A), where pf cs {A) is 
the BCS mean-field superfluid density with a modified 
pairing gap A. The effective gap is now renormalized 



to A 



(A s 



+ a 2 W 2 



where A pg is a temperature- 



dependent pseudogap describing the effect of bound pairs 
of the Fermi excitations. As shown in Fig. [TJ the pseu- 
dogap p s is very similar to the rescaled BCS result at 
higher temperatures. At low temperatures, the prefac- 
tor A 2 C /A 2 in the pseudogap expression for p s leads to 
a normal fluid density p n = p — p s oc T 3 / 2 (lik e an ideal 
Bose gas of molecules). We refer to Ref. [34j for more 
details. 

It is still not clear how to assess the treatment of 
bosonic pairing fluctuations used in the pseudogap cal- 
culation I33l 13411 . One indication of what it misses is to 



consider the BEC limit of the crossover, in which case 
the superfluid density predicted by Ref. [33| reduces to 
the condensate fraction of a noninteracting Bose gas of 
molecules, rather than the superfluid density for a gas of 
Bogoliubov excitations. This suggests that the pseudo- 
gap theory does not have any self-consistency problem 
near T c because it leaves out the interactions between 
bound pairs (equivalent to working with an ideal Bose 
gas). In future work, we will give a more detailed com- 
parison between the NSR-type fluctuation theory we use 
and the renormalized mean-field BCS theory involving a 
pseudogap. 



VI. BREATHING MODE FREQUENCIES 

The variational equations for the breathing mode fre- 
quencies using the ansatz in Eq. (|4T|) are @ 



ij w nj 



MniUJ 



ji/Q- n j + 2kj i a s j 



Here the "mass moments" Mi are defined by 

M si = J dr p sQ x 2 , M ni ee J dr p n0 x 2 , 
and the "spring constants" kfj , , and fcf™ are 



k n — 



dr 



dr 



dp \ djpsQXj) d(p s0 Xj) 
dp J s dxi dxj 

dp \ d(p n oXi) d(p n oXj) 

dp, 



(52) 



(53) 



(54) 



(55) 



dxi 



dx,i 



+ 



,'dT\ djpnQXj) d(s Xj} 
dp J dxi dxj 



dT \ d(s Xj) d(s Xj) 
ds J dxi dxj 



(56) 



and 



dr 



dp \ d(p S QXi) d(p n0 Xj) 
dp, 



dxi 



dx* 



dT \ d(p s0 Xi) d(s Xj) 
dp 



dxi 



dx-j 



(57) 



To solve for the breathing modes (see Section fVIip . it 
is useful to rewrite the above equations. We define the 
following new coefficients involving the spring constants: 



K^ = 2kt J+ 2k™ 



(58) 
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and 



(59) 



Adding the two equations for the breathing modes in 
Eqs. and (53]), we obtain 

uj 2 (M st a sl + M m a m ) = - ^ (K^a nj + K^a sj ) . 

j 

(60) 

Furthermore, dividing Eqs. ([52]) and ([53]) by M S i and 
M n i respectively, and subtracting one from the other, we 
obtain 



uj 2 (a si - a m ) = - 

j 



M sl 

K s 
M, 



2k( 

Mr, 



K n 



K s 

3 

M n 



I IV) 



(61) 



After some rearranging, we can write the coefficients 
defined in Eqs. (55]) and as 



K s = 2 



and 



dr 



dxi 



KV 



dr 



dr 



dspXj 
dx. 



Vp 



J Vxi 



Vp 

' T>Xj 



VT 



Po 



dp, 
dp' 



-•so 



+Po 



dp\ 



-so 



dp\ 

9s ) p 
(62) 



dp 



Vxj 



'dT\ fdT\ 



p 

(63) 



Here we have defined [not to be confused with the La- 
grangian derivative defined in Eq. (|16p ] 



Vp 


[dp" 


) ^ + 




Vxj " 




' s 9x 3 


m 



8x4 



(64) 



and 



VT _ 


fdT x 


) ^ + 




Vxj ~ 


\9 P/ 




m 



£ ^+ £ ^- (65) 



In writing down these equations, we also have made use 
of the Maxwell relation given by Eq. ([30]) . 

We next proceed to show that the expressions given 
in Eqs. (|62p and (|63p can be written in terms p s p, p n o, 
and the two thermodynamic derivatives, (dP/dp) Sl and 
(dP/ds)p. To handle the derivatives V(p,T)/Vxj, we 
note that the gradient of the equilibrium chemical po- 
tential po can be written as [using Eq. ([TO]) ] 



d PJs P ° \ ds 



V.s + W C5 



(66) 



Recall that in equilibrium, both the temperature and the 
chemical potential are spatially uniform (V/io = VT = 
0). Thus, for a harmonic trapping potential given by 
Eq. ©, Eq. §jo§ reduces to 



[dp' 


) ^ + 


m 


\dp, 







p UUj 3 



dx 



= -upj, (67) 



where u>j is the trap frequency along the rrj-axis. Us- 
ing the results in Eqs. (3J]) and (57]), Eqs. JM]) and (163)) 
simplify to 



Vp 

Vxj 



-^) X 3 



and 



VT 

Vxj 



= 0. 



(68) 



(69) 



Using Eqs. (3T]), (32]), (68]), and (69]), and integrat- 
ing by parts, the new spring constants KL and in 
Eqs. (S2) and (53]) reduce to 



2 dr p s0 Xi 



9A 2 d fdP 

28^ Xi - — 



and 



K n 



dr^p n0 x 
d fdP 



2S lj uj l Xi 



d 

dxi 



-SpXi 



dxi \ds J 



).}■ 



dp 



dP 
~dp~ 



(70) 



(71) 



In summary, we have reduced the algebraic equations 
for the breathing modes given by Eqs. ([52"]) and ([53]) to 
the set of equations given by Eqs. ([60]) and (|61[) . with 
the simpler spring constants given by Eqs. (TO"]) and (7TJ). 
Unlike the original spring constants defined in Eqs. (|55|) - 
(|57p . these new spring constants K^ n only involve deriva- 
tives of the pressure. 

Making use of the special properties of universal ther- 
modynamics at unitarity, the spring constants K?- and 
given in Eqs. (TDJ) and (7TJ) reduce to simple ex- 
pressions that involve only the mass moments M S i, M n i 
and the trap frequencies w,. Using Eqs. © and (ITS]) in 
Eqs. (70]) and (7T]), we obtain: 



Kf 



2M si (2%+2/3)o;? 



and 



K% =2M m (2%+2/3K 2 



(72) 



(73) 



These results (valid at unitarity) will be used in the next 
section. 
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VII. BREATHING MODES AT UNITARITY 



Using Eqs. ([72]) and ([73]) in Eq. (|60j( and (JHTJ), the 

variational equations reduce to 



uj 2 (M si a si + M„,a M 

J2( 26 a + 2 / 3 H 2 (M aJ a Bi + :u ni <i 



-■ , : ' n j" nj I ■ (74) 



and 



u) 2 (a s 



E 



Qui 
M,, 



M„ 



Ms. 
M„ 



(2*«+2/3)w? 



(75) 



These equations for the variational parameters a s , a n will 
be used to determine the hydrodynamic breathing modes 
at unitarity. 

By inspection, one immediately sees that Eq. ([75]) has 
a solution given by 

a S i = a ni . (76) 

This solution corresponds to a solution of the Landau 
two-fluid equations of the form 



v s (r,t) = v„(r,f). 



This describes the expected locally isentropic (also 
isothermal) breathing mode at unitarity. In Sec- 
tion IVII A[ we show that the frequency of this mode 
is independent of temperature, as argued by Thomas et 
al. 0| . Substituting the in-phase solution Eq. (f76"j) into 
Eq. ([74"]) , the latter reduces to 



MiLu 2 ai = ^2 ( 2S v + 2 / 3 ) a 



(78) 



where Mi = M s i + M n i and a S i — a n i = a,i. From the 
definition of the spring constants Kfj and if™- in Eqs. ([55]) 
and (59]), one sees that + K% = Kfc + Applying 
this result, we see that the expressions in Eqs. (172"]) and 
(1731) imply (valid at unitarity) 



MiUj 2 = MjLO 2 



(79) 



for all coordinates, i,j = x,y, z. Making use of this result 
in the right-hand side of Eq. ([75]) , it reduces to 



uj 2 1B a t = 2oj 2 a, + -uj 2 ^ a i> 

j 



(80) 



where uj\b denotes the frequency of the in-phase breath- 
ing mode. Equations (|76[) and ([80]) describe the in-phase 
oscillation of the normal and superfluid components cor- 
responding to v s — v„. We discuss the solutions of 
Eq. (JEHD below. 



In addition to the in-phase solution in Eq. ([76]) . there 
is an out-of-phase solution corresponding to 



M si a si + M ni a ni = 0, 



(81) 



which satisfies Eq. ([71]) . Substituting this out-of-phase 
solution into Eq. ([75"]) , we find a closed equation for the 
a s i parameters and the frequency ujib of the out-of-phase 
breathing mode at unitarity, namely 



J 2B Ll si 



_^ 2% + 2/3 )w 2 

M rj M rj M ni 1 ' 



(82) 

Here we have defined the reduced mass moment M r i as 
M ri = . (83) 



A. In-phase mode at unitarity 

The in-phase mode given by Eq. I|80p is a normal mode 
of the Landau two-fluid equations at unitarity, valid at 
all temperatures. It is equivalent to Eq. (3) in Ref. [35 1 
for the zero-temperature breathing mode frequency of a 
trapped Fermi gas at unitarity, assuming a polytropic 



(77) equation of state [22j /i(p) = p 1 with polytropic expo- 



nent 7 = 2/3. For an axisymmetric trap, uj x = uj y = uj±, 
the axial and longitudinal breathing modes are charac- 
terized by solutions of the form a x = a v . In this case, 
the solution of Eq. ([8H]) is well known [35j , 



J 1B 



+ ± ^(10u;i-8u;f) 2 +32u;iu;i. (84) 



We conclude that the frequency of the in-phase two- 
fluid hydrodynamic breathing mode at unitarity is in- 
dependent of temperature and equal to the zero temper- 
ature value. We further note that for an isotropic trap 
(uj x = u>y = u> z = luq), we have a x = a y = a Zl and the 
in-phase mode frequency in Eq. (|84p reduces to 



u\b = 2cl>o. 



(85) 



Castin [36[ has argued that there is an exact eigenstate of 
the isotropic trap Hamiltonian such that all atoms move 
with velocity v(r, t) = arcoscut at all temperatures, giv- 
ing rise to a temperature independent mode with fre- 
quency to — 2uj$. It is reassuring that two- fluid hydro- 
dynamics gives a result in agreement with this predic- 
tion [3?J ■ The temperature independence of the in-phase 
breathing mode is also consistent with the results of the 
direct numerical solution of the Landau two-fluid equa- 
tions reported in Ref. 0]. 

The fact that the in-phase breathing mode frequen- 
cies are independent of temperature is a consequence of 
the special thermodynamic properties at unitarity, and 
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is not expected to hold away from unitarity. In typ- 
ical experiments where the trap is highly anisotropic 
(u) z <C w±), the radial hydrodynamic breathing mode 
frequency [given by the upper branch of Eq. (|84j) ] is well- 
approximated by 



Wis 



(86) 



The predicted temperature independence of the in-phase 
breathing mode frequency is consistent with the experi- 
mental results of Thomas and coworkers @ . They found 
only a small difference (a few percent) between the mea- 
sured radial breathing mode frequency and Eq. (|86[) over 
a large temperature range, including well into the normal 
phase. 



B. Out-of-phase mode at unitarity 

We now discuss the out-of-phase breathing mode at 
unitarity. In the special limit of an isotropic trap, a S i — > 
a s and M sl -> Mf /3, M m -> M*/3, where 



M. 



a _ 



dr p s0 (r)r 2 



(87) 



= / drp n0 (r)r 2 



and 



Using these identities, Eq. 

k B 



simplifies to 



J 2B 



Mi 3 Mi 3 



(89) 



where we have defined the breathing mode spring con- 
stant, 



hi 



dr 



Oft 
dp 



[V • (rp s0 (v))Y 



(90) 



The reduced mass moment is now given by M.f = 
M?M%/{M?+M%). As follows from Eq. flST]), this out- 
of-phase mode corresponds to the following eigenvector: 



Mfa s + M^a„ = 0. 



(91) 



We now present numerical results for the out-of-phase 
breathing mode. From Eq. (|89|) . we see that only two 
thermodynamic functions enter in the evaluation of the 
mode frequency: the superfluid density p s (r) and the 
adiabatic compressibility (dp/dp) s (r). The calculation 
of these quantities in a uniform gas within an NSR-type 
formalism is discussed in Section [V] We use a local den- 
sity approximation (LDA) to calculate the local super- 
fluid density and compressibility in a trapped Fermi gas. 



The local density approximation in an isotropic har- 
monic trap (loq) amounts to determining the global chem- 
ical potential p from the local equilibrium condition 



P = Mhom \p(r),T/T F (p)] + uj 2 r 2 /2. 



(92) 



Here the local reduced temperature T/Tp{p) depends on 
the local mass density p(r) . Eq. (|9"2"j) is solved for the the 
density profile p{r), subject to the constraint 



drp(r) = Nm. 



(93) 



To solve for p(r) using LDA, for a given temperature T, 
we tabulate the local chemical potential as a function of 
the mass density using Eq. (pfB")) . With an initial guess 
of the global chemical potential, we determine the local 
chemical potential from the local equilibrium condition 
in Eq. (|92p . and invert it in tabular form to find the mass 
density. The global chemical potential is then adjusted 
slightly to enforce the number conservation requirement 
in Eq. (|93|) . giving a better estimate for the next iterative 
step. 

In a harmonic trap, it is convenient to use the trap 
units, where m = fcs = h = luq = 1, i.e., we take the 
characteristic harmonic oscillator length dho = \fh/mujQ 
and the characteristic level spacing Huiq as the units of 
the length and energy, respectively. We use the Fermi 
energy Ep = (3A r ) 1 / 3 fiw and the corresponding temper- 
ature Tp = Ep/ks of an ideal Fermi gas to characterize 
the energy scale and the temperature scale, where N is 
the total number of atoms. The distance and the mass 
density are conveniently given in units of the Thomas- 
Fermi radius Rtf — (24JV) 1 ^ 6 ai lo for an ideal Fermi 
gas, and the mass density at the centre of the trap, 
Ptf = (24N) 1 / 2 /(3TT 2 )ma^, respectively. 

In Fig. [2] we plot the profiles for the total mass density, 
the superfluid mass densities using the results in Fig. [TJ 
and the adiabatic compressibility. The total mass den- 
sity shows a bi-modal distribution, as expected from the 
general universal argument [2l[ ■ The superfluid densities 
drop to zero steeply at the superfluid-normal interface in 
the trap. 

Having calculated the local adiabatic compressibility 
and the superfluid and total mass density profiles, it is 
straightforward to evaluate the mass moments and the 
spring constant that enter the out-of-phase breathing 
mode frequency. 

The frequency ll>2b is plotted in Fig. [3] using two ap- 
proximations for the superfluid density, as a function of 
temperature. One immediately sees that the out-of-phase 
breathing mode frequency is quite sensitive to the su- 
perfluid density in a uniform gas (using the LDA). This 
underlines the importance of calculating p s with better 
accuracy as input into our variational theory. However, 
the qualitative features of the temperature dependence 
of oj 2 b ar e similar for both the fitted NSR and scaled 
mean- field BCS data for p s . Namely, the frequency in- 
creases rapidly at low temperatures and decreases with 
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FIG. 2: (color online) Profiles of several thermodynamic func- 
tions at temperature T = 0.20Tf, where TV is the Fermi 
temperature of a trapped ideal Fermi gas. Here we have 
used the results given by the NSR-type Gaussian fluctua- 
tion theory The superfluid transition temperature is 
T c ~ 0.27Tf. The adiabatic compressibility is only needed in 
the superfluid region of the trap. 
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FIG. 3: The out-of-phase breathing mode at unitarity as a 
function of temperature for an isotropic trap using an im- 
proved NSR calculation for the equation of state [123 ] ■ The 
frequencies obtained using the fitted NSR p B data and the 
scaled BCS mean-field data are given by the solid and dashed 
lines, respectively (see Fig. (TJ. The arrow indicates the su- 
perfluid transition temperature, T c ~ 0.27Tf. 



In Appendix [Bl we argue that the divergence of u>2B as 
T — ► is not an artifact of the local density approxima- 
tion (LDA) we use to evaluate the coefficients in Eq. 
We compare the results of a mean-field LDA calcula- 
tion directly with the results obtained by self-consistently 
solving the Bogoliubov-de Gennes equations and find ex- 
cellent agreement. The discussion in Appendix [Bl shows 
that LDA is not the source of any significant error at low 
temperatures. 

The increase of loib as T — > can be understood within 
our variational formalism as follows. As emphasized in 
Ref. [5[, our variational solutions of the two-fluid equa- 
tions describe two coupled harmonic oscillators with ef- 
fective masses given by the mass moments for the mode 
in question [for the breathing mode, these are given by 
Eqs. (J57]) and (EH])]- As T — ► 0, the mass of the normal 
fluid "oscillator" goes to zero. As with two coupled har- 
monic oscillators, in this limit, the small (normal fluid) 
mass executes a high frequency (and large amplitude) 
oscillation about the heavy (superfluid) mass, which is 
essentially static. We should also caution that at low but 
finite T, the Landau two-fluid equations are no longer 
valid because local equilibrium cannot be established. 

The high-temperature (T — > T c ) behavior of u>2B is 
discussed in AppendixObased on the analytic expression 
in the BCS approximation for the superfluid density in a 
trapped gas [381 ]. 



VIII. DIPOLE MODES 

The dipole modes discussed in Ref. ||| are character- 
ized by the uniform displacement fields given by Eq. ([43[) . 
Inserting this ansatz into the action given in Eq. (f38|) 
and taking its variation, one finds an in-phase oscillation 
(generalized Kohn mode) with a s = a n and frequency 
u>id = w z given by the trap frequency lo z along the z-axis. 
In addition, there is an out-of-phase mode corresponding 
to the solution M s a s + M n a n = 0. The frequency of this 
mode is given by 0] 



J 2D 



MP ' 



(94) 



Here M D = M D M D / ( M D + M D ) is the reduced mass 
of the superfluid and normal fluid components, with 



Mi 



dr p s0 (r) 



(95) 



increasing temperature, before increasing again as T c is 
approached. In both cases, the frequency of the out-of- 
phase breathing mode is larger than the in-phase breath- 
ing mode uiib = 2w - These features are quite differ- 
ent from the results of He et al. Q- They find the out- 
of-phase breathing mode frequency starts below the in- 
phase mode frequency at low temperatures, and increases 
monotonically as T approaches T c . 



and 



M. 



D 



dr p n o(r) 



(96) 



giving the masses of the superfluid and normal fluids. 
The spring constant k® n in Eq. (f94|) is defined as 0] 



, D = 
v sn — 



dr 



dp J s dz 



dT\ dso 
dp) s dz 



dz ' 



(97) 
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FIG. 4: The frequency of the out-of-phase dipole mode at 
unitarity in an isotropic trap as a function of temperature. 
See caption of Fig. |3] 



This is the analogue of the corresponding spring constant 
kfj 1 for the breathing mode, defined in Eq. ([57j) . We note 
that k® n is the negative of the analogous spring constant 
k sn in Eq. (79) of Ref. [|. 

It is convenient to write the out-of-phase mode fre- 
quency in Eq. (|94[) in terms of a simpler spring constant 

, which only involves the isentropic compressibility 
(dp/dp) s . Using Eqs. ^30j) and (J67|) in Eq. (J97]), we find 



- / dr 



where 



k? = I dr 



Oft 
dp 



dpso 



dp s o 
dz 



(98) 



dp 
dp 



( dpso 
V dz 



(99) 



Using these results in Eq. (|94| . we find it reduces to 



J 2D 



Ms 
M, 



(100) 



This frequency of the out-of-phase dipole mode in an 
isotropic trap (u> z = u>q) is plotted in Fig. 0] as a func- 
tion of temperature. As in Fig. [3], we compare the results 
obtained using the superfluid density given by NSR the- 
ory and the BCS mean-field approximation (see Fig. [!}. 
We note that the expression given in Eq. (|100|) for the 
frequency of the out-of-phase dipole mode is very simi- 
lar to the formula for the frequency of the out-of-phase 
isotropic breathing mode at unitarity given in Eq. (|89[> . 
Thus it is not surprising that the frequencies (shown in 
Figs. [3] and 2]) of both of these out-of-phase modes exhibit 
similar behavior as a function of temperature. 



As discussed in Section HVi the in-phase dipole mode 
considered in this Section is also an example of a locally 
isentropic mode. Unlike the in-phase breathing mode, 
however, which is only isentropic at unitarity, the gener- 
alized Kohn mode is always characterized by v s = v„ , at 
all temperatures and interaction strengths. This feature 
also follows from the condition given in Eq. Q35p since 
V • u = for the dipole mode. 



IX. CONCLUDING REMARKS 

In this paper, we have presented results for the breath- 
ing and dipole mode solutions of the Landau two-fluid 
equations for a Fermi superfluid at unitarity in an 
isotropic trap. Our work is based on a recent varia- 
tional formulation [5[ of the two-fluid equations. We have 
shown that the variational equations simplify at unitar- 
ity, where the coefficients only depend on the compress- 
ibility and the superfluid density. Understanding the na- 
ture of Fermi gases at unitarity, where the s-wave scatter- 
ing length diverges, is a challenging many-body problem. 
In contrast to the in-phase dipole and breathing modes, 
which we have shown to be independent of temperature, 
the out-of-phase modes are very dependent on tempera- 
ture. Measurement of the out-of-phase mode frequencies 
will provide a sensitive test of current microscopic theo- 
ries of a Fermi gas at unitarity, including the predictions 
of "universal thermodynamics" pl| . In particular, our 
mode frequencies in Figs. [3] and 2] show that the results 
are very dependent on the temperature dependence of the 
superfluid density. It would be very useful to have a more 
accurate ab-initio calculation of p s (such as Ref. [30]|). 

In a companion paper [To| . we show that the fre- 
quencies of these hydrodynamic modes can be measured 
using two-photon Bragg spectroscopy, a standard tool 
used to study excitations in trapped, ultracold quantum 
gases [TTj ] . 

We emphasize that the results presented in this paper 
for the frequencies of the low-lying out-of-phase breath- 
ing and dipole modes are based on the simplest possible 
variational ansatz for these modes. Our variational re- 
sults provide an upper bound on the exact frequencies 
for these modes. In future work, we will discuss results 
based on an improved variational ansatz. 

While we have concentrated on the two-fluid modes at 
unitarity in the present paper, our general variational 
formulation of the two-fluid equations can be applied 
anywhere in the BCS-BEC crossover, as long as the in- 
teractions are sufficiently strong to ensure collisionally 
hydrodynamic behavior. It would be interesting to con- 
sider the two-fluid modes of a strongly-interacting Bosc- 
condensate of dimer molecules, on the BEC side of uni- 
tarity. 

We can use the thermodynamic functions discussed in 
Section[V]to evaluate the temperature dependence of first 
and second sound velocities at unitarity in a uniform gas 
based on the NSR theory. However, the physics of second 
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sound in a uniform gas [20| is quite different from the out- 
of-phase modes in a trapped gas discussed in the present 
paper. One finds that at unitarity, the BCS Fermi ex- 
citations make the dominant contribution to the second 
sound velocity but as one goes to lower temperature, the 
increasing gap Ao (see, for example, Ref. [l!|) freezes 
out this contribution relative to the undamped bosonic 
excitations. The end result is that as T — > 0, the second 
sound velocity increases and approaches c/v3) where c 
is the Bogoliubov phonon velocity. We will give a more 
complete discussion of second sound in a uniform gas in 
another publication. 
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APPENDIX A: FIRST SOUND IN SUPERFLUID 
4 HE AS A LOCALLY ISENTROPIC MODE 



(see Sec. 16 in Landau and Lifshitz [39j ]) 



In part A of Section IVII1 we showed that at unitarity 
the in-phase breathing mode is a locally isentropic mode 
(VST = 0), where the local superfluid and normal fluid 
velocities are equal, v s (r, t) = v„ (r, t) . It is useful to 
compare this analysis with the case of superfluid 4 He, 
where first sound also describes a locally isentropic mode. 

The two-fluid modes in uniform superfluid helium are 
to a very good approximation given by [l|, Q v s = v n 
and psov s + p n oVn = 0, corresponding to first and second 
sound, respectively. First sound describes a locally isen- 
tropic density oscillation (ST — 0), while second sound 
describes a pure temperature (5p — 0) oscillation 
Note that in a uniform superfluid, the condition 'VST = 
is equivalent to ST = since a uniform oscillation of the 
temperature is impossible. 

The existence of a locally isentropic first sound mode in 
uniform 4 He is also accounted for by the condition we give 
in Eq. (|35[) . In a uniform system, all equilibrium ther- 
modynamic quantities are independent of position and 
the term in the second line Eq. (|35[) that involves the 
gradient of (dP/ds) p vanishes (recall that in a trapped 
superfluid, it only vanishes at unitarity). Since Eq. (|3"T|) 
is not satisifed by the plane-wave solutions of the uni- 
form two-fluid equations, we see that the condition for 
a locally isentropic first sound mode to exist is given by 
Eq. (|3"6"|) , namely that (dP/ds) p = 0. Using the identity 



dP 
~d7 



T_ /dP\ 
PCv \dT) p 



(Al) 



where c v — T(ds/dT) p is the equilibrium specific heat 
per unit mass, one sees that (dP/ds) p ~ implies 
(dP/dT)p ~ 0. In this case, the adiabatic and isother- 
mal compressibilities are equal {(dP/dp) s — (dP/dp)r]- 
When dealing with the two-fluid equations in superfluid 
helium, this equivalence leads to a well known simpli- 
fication in the equations for first and second sound [2j . 
The simplified equations can be easily solved leading to 
the result that first sound is a locally isentropic mode 
(v s = v„), with a sound speed given by the adiabatic 
compressibility Ux = y/ (dP/dp)g. 



APPENDIX B: FREQUENCIES CLOSE TO T = 
IN THE BCS APPROXIMATION 



In this Appendix, we show that the local density ap- 
proximation is not responsible for the diverging frequency 
of the out-of-phase breathing and dipole modes as T — > 
(see Figs. [3] and E}. In this limit, the reduced mass mo- 
ment is approximately the mass moment of the normal 
component, M r — > M n . Thus, the out-of-phase mode 
frequency is inversely proportional to the normal mass 
moment, which becomes very small at low temperature. 
One may question the numerical accuracy of the calcula- 
tions. In particular, the strong temperature dependence 
of the out-of-phase mode frequencies at low T might be 
an artifact of the local density approximation used in 
Figs. [3] and 2] To check this point, we calculate the 
out-of-phase breathing mode frequency using the ther- 
modynamic functions for a trapped gas given by directly 
solving the Bogoliubov-de Gennes (BdG) equations for a 
finite number of atoms. 

We solve the coupled BdG equations for the Bogoli- 
ubov quasiparticles of a Fermi gas in an isotropic har- 
monic trap at unitarity. A microscopic expression of the 
superfluid density of a finite size inhomogeneous system 
may be derived by considering the moment of inertia of 
the Fermi gas, or equivalently, by calculating the increase 
in free energy after impo sing a twisted boundary phase 
for the order parameter [18j |. In an isotropic trap, the 
superfluid density is given by 



Pso (r) 



Po (r) 
1(1 + 1) (21 + 1) 



nl 



df(E nl 



8tt 



Ki (r) + v 2 nl (r)] . (Bl) 



Here u n i (r) and v n i (r) are the radial wavefunctions 
of the Bogoliubov quasiparticles. The full wavefunc- 
tions have the form Uj (r) = [u ni (r) /r] Y tm (6, <j>) and 
Vj (r) = [v nl (r) /r] Y lm (0, <f>). The quasiparticle energy 
E n i is allowed to be negative, and the summation of the 
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FIG. 5: Comparison of an LDA calculation of the frequency of 
the out-of-phase breathing mode with a calculation using the 
Bogoliubov-de Gennes equations. Here the LDA result (solid 
line) is obtained by using the BCS mean-field theory for a 
Fermi gas at unitarity. The solid squares are calculated us- 
ing the self-consistent Bogoliubov-de Gennes mean-field equa- 
tions. 



level indices (nl) is over both positive and negative energy 
levels. / (x) is the Fermi-Dirac distribution function. 

Fig. [5] compares the frequency of the out-of-phase 
breathing mode calculated (a) within LDA using a mean- 
field BCS equation of state and (b) from a self-consistent 
calculation of the BdG equations. The BdG frequency 
(for N = 2 x 10 5 atoms) is smaller but close to the LDA 
result. The increase of the mode frequency with decreas- 
ing temperature is clearly seen in both the BdG and LDA 
calculations. The disagreement seen in Fig. [5] is not un- 
expected since numerical calculations we carried out as a 
function of N show that the BdG results converge very 
slowly with increasing N. We conclude that the strong 
increase of the out-of-phase mode frequency in the low 
temperature regime is not an artifact of the LDA. 



der parameter in a harmonic trap [36 



A(r) cx T c 



T c -T 



1/2 



'-a 



1/2 



(CI) 



Here the radius R c cx ^J5T/T C and 5T = T c - T. Since 
the BCS supcrfluid density varies as p s o (r) cx A 2 (r) near 
T c , we may write 



p s0 (r) = a 1 - 



Rl 



(C2) 



where the prefactor a cx 8T/T C . We note that even 
though both a and R\ vanish linearly with temperature 
close to T c , the ratio ajR 2 c remains finite. 

In the vicinity of T c , the superfluid mass moment is 
much smaller than the normal mass moment, and the 
reduced mass moment of both the dipole and breathing 
modes reduces to the superfluid mass moment, M r — > 
M 8 . Taking the breathing mode as an example, its mode 
frequency [given by Eq. (|8"9")) ] reduces to U)\ B = kf /Mf 
(note that kf /Mf remains finite as T — ► T c , while 
Mf /Mf vanishes). Since R c <C 1, the adiabatic com- 
pressibility is nearly constant in the region of interest, 
and we denote it as 7 = (dp/dp) s for T — > T c . 

The calculations of the superfluid mass moment and 
the spring constant are straightforward. Substituting Eq. 
(fC2j) into Eqs. (87} and (90]), we obtain 



and 



16tt 



(C3) 



(C4) 



Thus, using a weak-coupling BCS mean-field calcula- 
tion near T c , the frequency of the out-of-phase breathing 
mode is predicted to be 



APPENDIX C: FREQUENCIES CLOSE TO T c IN 
THE BCS APPROXIMATION 

Close to the superfluid transition temperature, the 
temperature dependence of the out-of-phase mode fre- 
quencies using the mean-field BCS superfluid density 
can be worked out without using the LDA. In this re- 
gion, an analytic result for the weak-coupling BCS super- 
fluid density in a trapped gas is given by Baranov and 
Petrov [38| . Assuming a second order phase transition 
near T c , Ginzburg-Landau theory predicts the following 
temperature dependence for the position dependent or- 



J 2B 



MP 



(C5) 



As noted earlier, both a/R 2 and 7 approach constant val- 
ues close to T c , Thus, the out-of-phase breathing mode 
frequency is finite at the transition temperature. We have 
checked the validity of Eq. (|C5[) by numerically calculat- 
ing the values of a, 7, and R c . As these parameters do 
not change much above the temperature 0-5T c ^ ra p, the 
mode frequency becomes fairly constant in this temper- 
ature range, in agreement with our LDA results in Fig. [3] 
for the breathing-mode frequency based on the scaled 
BCS superfluid density (dashed curve). 
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